## Create simulated data many times

simmer <- function(meth="ktimes", kfac = 2, stat="mean", covs="cont", ncovsexact=0, ncovscont = 2, mu =0, sd = 1, r=0, restr.list = NULL, n.tr=2, whenoutlier = NULL, howoutlier = "10timesmax", units=40, times=100, outfileDir=""){
	
	for (idx.times in 1:times){
	
	## make bivariate normal data
	if(length(mu) < ncovscont){
		mmm <- rep(mu, ncovscont)
	}else{
		mmm <- mu
	}
	
	if(!(is.matrix(r))){
		cormat <- diag(1, ncovscont)
		cormat[row(cormat)!=col(cormat)] <- r
		SSS <- (sd^ncovscont)*cormat
		x <- mvrnorm(units, mmm, SSS)
	}else{
		xstan <- mvrnorm(units, rep(0, nrow(r)), r)
		x <- round(xstan*sd+mmm)
	}
	
	## draw exacts from unif(cats)
	if(ncovsexact>0){
		for(idx.exdraw in ncovsexact:1){
			tmp.exdraw <- sample(restr.list[[idx.exdraw]], units, rep=TRUE)
			x <- cbind(tmp.exdraw, x)
			colnames(x)[1] <- cov.names[idx.exdraw]
		}
	}
		
	## check draws given restrictions
	if(!(is.null(restr.list)) && (length(restr.list)>ncovsexact)){
		for(idx.unit in 1:units){
			for(idx.var in 1:ncol(x)){
				if(!(x[idx.unit, idx.var] %in% restr.list[[idx.var]])){
					x[idx.unit, idx.var] <- sample(restr.list[[idx.var]],1)
				}
			}
		}
	}

	if(covs=="contbim"){
		## make bimodal from normal data
		mmm <- rep(mu, ncovscont)
		cormat <- diag(1, ncovscont)
		cormat[row(cormat)!=col(cormat)] <- r
		SSS <- (sd^ncovscont)*cormat
    consider.factor <- 50
		x <- mvrnorm(units*consider.factor, mmm, SSS)
		x <- x[order(x[,1]),]
		x <- x[c(1:(units/2), (units*consider.factor-units/2+1):(units*consider.factor)),]
	}
	
	cvnames <- paste("cont", 1:ncovscont, sep="")
	
	if(ncovsexact>0){
		exact.v <- cov.names[ncovsexact]
		exact.values <- x[1,1:ncovsexact]
	}else{
		exact.v <- exact.values <- NULL
	}
	
	if(!is.null(whenoutlier)){
		if(whenoutlier == 1){
			if(howoutlier == "10timesmax"){
				seqblock1(id.vars = "id", id.vals =1, exact.vars = exact.v, exact.vals = exact.values, covar.vars = cvnames, covar.vals=x[1,(ncovsexact+1):ncol(x)]*10, n.tr = n.tr, assg.prob.stat = stat, assg.prob.method = meth, assg.prob.kfac = kfac, file.name = paste(outfileDir, "/sim",idx.times,".RData", sep=""), verbose=F)
			}
			if(howoutlier == "5fixed"){
				seqblock1(id.vars = "id", id.vals =1, exact.vars = exact.v, exact.vals = exact.values, covar.vars = cvnames, covar.vals=rep(5, ncol(x)-(ncovsexact+1)), n.tr = n.tr, assg.prob.stat = stat, assg.prob.method = meth, assg.prob.kfac = kfac, file.name = paste(outfileDir, "/sim",idx.times,".RData", sep=""), verbose=F)
			}

		}else{
			seqblock1(id.vars = "id", id.vals =1, exact.vars = exact.v, exact.vals = exact.values, covar.vars = cvnames, covar.vals=x[1,(ncovsexact+1):ncol(x)], n.tr = n.tr, assg.prob.stat = stat, assg.prob.method = meth, assg.prob.kfac = kfac, file.name = paste(outfileDir, "/sim",idx.times,".RData", sep=""), verbose=F)		}			
	}else{
		seqblock1(id.vars = "id", id.vals =1, exact.vars = exact.v, exact.vals = exact.values, covar.vars = cvnames, covar.vals=x[1,(ncovsexact+1):ncol(x)], n.tr = n.tr, assg.prob.stat = stat, assg.prob.method = meth, assg.prob.kfac = kfac, file.name = paste(outfileDir, "/sim",idx.times,".RData", sep=""), verbose=F)
	}
	if(!is.null(whenoutlier)){
		if(whenoutlier > 2){
			for(idx.units in 2:(whenoutlier-1)){
				seqblock2k(object= paste(outfileDir, "/sim", idx.times, ".RData", sep=""), id.vals=idx.units, covar.vals = x[idx.units,], file.name = paste(outfileDir, "/sim", idx.times, ".RData", sep=""), verbose=F)			}
		}	
		for(idx.units in whenoutlier){
			if(howoutlier == "10timesmax"){
				seqblock2k(object= paste(outfileDir, "/sim", idx.times, ".RData", sep=""), id.vals=idx.units, covar.vals = x[idx.units,]*10, file.name = paste(outfileDir, "/sim", idx.times, ".RData", sep=""), verbose=F)
			}
			if(howoutlier == "5fixed"){
				seqblock2k(object= paste(outfileDir, "/sim", idx.times, ".RData", sep=""), id.vals=idx.units, covar.vals = c(5,5), file.name = paste(outfileDir, "/sim", idx.times, ".RData", sep=""), verbose=F)
			}						
		}
		for(idx.units in (whenoutlier+1):units){
			seqblock2k(object= paste(outfileDir, "/sim", idx.times, ".RData", sep=""), id.vals=idx.units, covar.vals = x[idx.units,], file.name = paste(outfileDir, "/sim", idx.times, ".RData", sep=""), verbose=F)
		}
	}else{for(idx.units in 2:units){
		seqblock2k(object= paste(outfileDir, "/sim", idx.times, ".RData", sep=""), id.vals=idx.units, covar.vals = x[idx.units,], file.name = paste(outfileDir, "/sim", idx.times, ".RData", sep=""), verbose=F)}
	}
  }	
}
